home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / bsp_knot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  17.2 KB  |  523 lines
  1. /******************************************************************************
  2. * Bsp_knot.c - Various bspline routines to handle knot vectors.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. #define KNOT_IS_THE_SAME 1e-10
  13.  
  14. /******************************************************************************
  15. * Returns TRUE iff The given knot vector has open end conditions.          *
  16. ******************************************************************************/
  17. CagdBType BspKnotHasOpenEC(CagdRType *KnotVector, int Len, int Order)
  18. {
  19.     int i,
  20.         LastIndex = Len + Order - 1;
  21.  
  22.     for (i = 1; i < Order; i++)
  23.        if (!APX_EQ(KnotVector[0], KnotVector[i]))
  24.         return FALSE;
  25.  
  26.     for (i = LastIndex - 1; i >= Len; i--)
  27.        if (!APX_EQ(KnotVector[LastIndex], KnotVector[i]))
  28.         return FALSE;
  29.  
  30.     return TRUE;
  31. }
  32.  
  33. /******************************************************************************
  34. * Returns TRUE iff t is in the parametric domain as define by the knot vector *
  35. * KnotVector its length Len and the order Order.                  *
  36. ******************************************************************************/
  37. CagdBType BspKnotParamInDomain(CagdRType *KnotVector, int Len, int Order,
  38.                                    CagdRType t)
  39. {
  40.     CagdRType
  41.     r1 = KnotVector[Order - 1],
  42.     r2 = KnotVector[Len];
  43.  
  44.     return (r1 < t || APX_EQ(r1, t)) && (r2 > t || APX_EQ(r2, t));
  45. }
  46.  
  47. /******************************************************************************
  48. * Return the index of the last knot which is less than or equal to t in the   *
  49. * given knot vector KnotVector of length Len. APX_EQ is used for equality.    *
  50. * Parameter t is assumed to be in the parametric domain for the knot vector.  *
  51. ******************************************************************************/
  52. int BspKnotLastIndexLE(CagdRType *KnotVector, int Len, CagdRType t)
  53. {
  54.     int i;
  55.  
  56.     for (i = 0;
  57.      i < Len && (KnotVector[i] <= t || APX_EQ(KnotVector[i], t));
  58.      i++);
  59.  
  60.     return i - 1;
  61. }
  62.  
  63. /******************************************************************************
  64. * Return the index of the last knot which is less than t in the given knot    *
  65. * vector KnotVector of length Len.                          *
  66. * Parameter t is assumed to be in the parametric domain for the knot vector.  *
  67. ******************************************************************************/
  68. int BspKnotLastIndexL(CagdRType *KnotVector, int Len, CagdRType t)
  69. {
  70.     int i;
  71.  
  72.     for (i = 0;
  73.      i < Len && KnotVector[i] < t && !APX_EQ(KnotVector[i], t);
  74.      i++);
  75.  
  76.     return i - 1;
  77. }
  78.  
  79. /******************************************************************************
  80. * Return the index of the first knot which is greater than t in given knot    *
  81. * vector KnotVector of length Len.                          *
  82. * Parameter t is assumed to be in the parametric domain for the knot vector.  *
  83. ******************************************************************************/
  84. int BspKnotFirstIndexG(CagdRType *KnotVector, int Len, CagdRType t)
  85. {
  86.     int i;
  87.  
  88.     for (i = Len - 1;
  89.      i >= 0 && KnotVector[i] > t && !APX_EQ(KnotVector[i], t);
  90.      i--);
  91.  
  92.     return i + 1;
  93. }
  94.  
  95. /******************************************************************************
  96. * Return a uniform float knot vector for Len Control points and order Order.  *
  97. * If KnotVector is NULL it is being allocated dynamically.              *
  98. ******************************************************************************/
  99. CagdRType *BspKnotUniformFloat(int Len, int Order, CagdRType *KnotVector)
  100. {
  101.     int i;
  102.     CagdRType *KV;
  103.  
  104.     if (KnotVector == NULL)
  105.     KnotVector = KV = (CagdRType *) IritMalloc(sizeof(CagdRType) *
  106.                                 (Len + Order));
  107.     else
  108.     KV = KnotVector;
  109.  
  110.     for (i = 0; i < Len + Order; i++)
  111.     *KV++ = (CagdRType) i;
  112.  
  113.     return KnotVector;
  114. }
  115.  
  116. /******************************************************************************
  117. * Return a uniform open knot vector for Len Control points and order Order.   *
  118. * If KnotVector is NULL it is being allocated dynamically.              *
  119. ******************************************************************************/
  120. CagdRType *BspKnotUniformOpen(int Len, int Order, CagdRType *KnotVector)
  121. {
  122.     int i, j;
  123.     CagdRType *KV;
  124.  
  125.     if (KnotVector == NULL)
  126.     KnotVector = KV = (CagdRType *) IritMalloc(sizeof(CagdRType) *
  127.                                 (Len + Order));
  128.     else
  129.     KV = KnotVector;
  130.  
  131.     for (i = 0; i < Order; i++)
  132.     *KV++ = 0.0;
  133.     for (i = 1, j = Len - Order; i <= j; )
  134.     *KV++ = (CagdRType) i++;
  135.     for (j = 0; j < Order; j++)
  136.     *KV++ = (CagdRType) i;
  137.  
  138.     return KnotVector;
  139. }
  140.  
  141. /******************************************************************************
  142. * Apply an affine transformation to the given knot vector. Note affine        *
  143. * transformation for the knot vector does not change the spline.          *
  144. * Knot vector is translated by Translate amount and scaled by Scale.
  145. ******************************************************************************/
  146. void BspKnotAffineTrans(CagdRType *KnotVector, int Len,
  147.                     CagdRType Translate, CagdRType Scale)
  148. {
  149.     int i;
  150.     CagdRType KV0 = KnotVector[0];
  151.  
  152.     for (i = 0; i < Len; i++)
  153.     KnotVector[i] = (KnotVector[i] - KV0) * Scale + KV0 + Translate;
  154. }
  155.  
  156. /******************************************************************************
  157. * Creates an identical copy of a given knot vector KnotVector of length Len.  *
  158. ******************************************************************************/
  159. CagdRType *BspKnotCopy(CagdRType *KnotVector, int Len)
  160. {
  161.     CagdRType
  162.     *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len);
  163.  
  164.     CAGD_GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * Len);
  165.  
  166.     return NewKnotVector;
  167. }
  168.  
  169. /******************************************************************************
  170. * Reverse a knot vector of length Len.                          *
  171. ******************************************************************************/
  172. CagdRType *BspKnotReverse(CagdRType *KnotVector, int Len)
  173. {
  174.     int i;
  175.     CagdRType
  176.     *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len),
  177.     t = KnotVector[Len - 1] + KnotVector[0];
  178.  
  179.     for (i = 0; i < Len; i++)
  180.     NewKnotVector[i] = t - KnotVector[Len - i - 1];
  181.  
  182.     return NewKnotVector;
  183. }
  184.  
  185. /******************************************************************************
  186. * Returns all the knots in KnotVector1 not in KnotVector2.              *
  187. * NewLen is set to new KnotVector length.                      *
  188. ******************************************************************************/
  189. CagdRType *BspKnotSubtrTwo(CagdRType *KnotVector1, int Len1,
  190.                CagdRType *KnotVector2, int Len2,
  191.                int *NewLen)
  192. {
  193.     int i = 0,
  194.     j = 0;
  195.     CagdRType
  196.     *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len1),
  197.     *t = NewKnotVector;
  198.  
  199.     *NewLen = 0;
  200.     while (i < Len1 && j < Len2) {
  201.     if (APX_EQ(KnotVector1[i], KnotVector2[j])) {
  202.         i++;
  203.         j++;
  204.     }
  205.     else if (KnotVector1[i] > KnotVector2[j]) {
  206.         j++;
  207.     }
  208.     else {
  209.         /* Current KnotVector1 value is less than current KnotVector2: */
  210.         *t++ = KnotVector1[i++];
  211.         (*NewLen)++;
  212.     }
  213.     }
  214.  
  215.     return NewKnotVector;
  216. }
  217.  
  218. /******************************************************************************
  219. * Merge two knot vector KnotVector1 and KnotVector2 of length Len1 and Len2   *
  220. * respectively into one. If Mult is not zero then knot multiplicity is tested *
  221. * not to be bigger than Mult value. NewLen is set to new KnotVector length.   *
  222. ******************************************************************************/
  223. CagdRType *BspKnotMergeTwo(CagdRType *KnotVector1, int Len1,
  224.                CagdRType *KnotVector2, int Len2,
  225.                int Mult, int *NewLen)
  226. {
  227.     int i = 0,
  228.     j = 0,
  229.     k = 0,
  230.         m = 0,
  231.     Len = Len1 + Len2;
  232.     CagdRType t,
  233.     *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len);
  234.  
  235.     if (Mult == 0)
  236.     Mult = Len + 1;
  237.  
  238.     while (i < Len1 && j < Len2) {
  239.     if (KnotVector1[i] < KnotVector2[j]) {
  240.         /* Use value from KnotVector1: */
  241.         t = KnotVector1[i++];
  242.     }
  243.     else {
  244.         /* Use value from KnotVector2: */
  245.         t = KnotVector2[j++];
  246.     }
  247.  
  248.     if (k == 0 || (k > 0 && !APX_EQ(NewKnotVector[k - 1], t))) {
  249.         NewKnotVector[k++] = t;
  250.         m = 1;
  251.     }
  252.     else if (m < Mult) {
  253.         NewKnotVector[k++] = t;
  254.         m++;
  255.     }
  256.     }
  257.  
  258.     while (i < Len1)
  259.     NewKnotVector[k++] = KnotVector1[i++];
  260.     while (j < Len2)
  261.     NewKnotVector[k++] = KnotVector1[j++];
  262.  
  263.     /* It should be noted that k <= Len so there is a chance some of the new */
  264.     /* KnotVector space will not be used (it is not memory leak!). If you    */
  265.     /* really care that much - copy it to a new allocated vector of size k   */
  266.     /* and return it while freeing the original of size Len.             */
  267.     *NewLen = k;
  268.  
  269.     return NewKnotVector;
  270. }
  271.  
  272. /******************************************************************************
  273. * Creates a new vector with the given KnotVector Node values.              *
  274. * The nodes are the approximated parametric value associated with the each    *
  275. * control point. Therefore for a knot vector of length Len and order Order    *
  276. * there are Len - Order control points and therefore nodes.              *
  277. * Nodes are defined as (k = Order - 1 or degree):                  *
  278. *                                              *
  279. *      i+k                                      *
  280. *       -                 First Node N(i = 0)              *
  281. *       \                                      *
  282. *          / KnotVector(j)        Last Node  N(i = Len - k - 2)          *
  283. *       -                                      *
  284. *        j=i+1                                      *
  285. * N(i) = -----------------                              *
  286. *            k                                  *
  287. ******************************************************************************/
  288. CagdRType *BspKnotNodes(CagdRType *KnotVector, int Len, int Order)
  289. {
  290.     int i, j,
  291.     k = Order - 1,
  292.     NodeLen = Len - Order;
  293.     CagdRType
  294.     *NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * NodeLen);
  295.  
  296.     for (i = 0; i < NodeLen; i++) {
  297.     NodeVector[i] = 0.0;
  298.     for (j = i + 1; j <= i + k; j++)
  299.         NodeVector[i] += KnotVector[j];
  300.     NodeVector[i] /= k;
  301.     }
  302.  
  303.     return NodeVector;
  304. }
  305.  
  306. /******************************************************************************
  307. * Creates a new vector with t inserted as a new knot. Attempt is made to make *
  308. * sure t is in the knot vector domain.                          *
  309. * No test is made for the current multiplicity for knot t in KnotVector.      *
  310. ******************************************************************************/
  311. CagdRType *BspKnotInsertOne(CagdRType *KnotVector, int Order, int Len,
  312.                                 CagdRType t)
  313. {
  314.     int Mult = BspKnotFindMult(KnotVector, Order, Len, t) + 1;
  315.  
  316.     return BspKnotInsertMult(KnotVector, Order, &Len, t, Mult);
  317. }
  318.  
  319. /******************************************************************************
  320. * Inserts Mult knots with value t into the knot vector KnotVector.          *
  321. * Attempt is made to make sure t in knot vector domain.                  *
  322. * If a knot equal to t (up to APX_EQ) already exists with multiplicity i only *
  323. * Mult - i knot are been inserted into the new knot vector.              *
  324. * Len is updated to the resulting knot vector.                      *
  325. * Note it is possible to DELETE a knot using this routine by specifying       *
  326. * multiplicity less then current multiplicity!                      *
  327. ******************************************************************************/
  328. CagdRType *BspKnotInsertMult(CagdRType *KnotVector, int Order, int *Len,
  329.                                 CagdRType t, int Mult)
  330. {
  331.     int i, CurrentMult, NewLen, FirstIndex;
  332.     CagdRType *NewKnotVector;
  333.  
  334.     if (!BspKnotParamInDomain(KnotVector, *Len, Order, t))
  335.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  336.  
  337.     CurrentMult = BspKnotFindMult(KnotVector, Order, *Len, t);
  338.     NewLen = *Len + Mult - CurrentMult;
  339.     NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) *
  340.                                 (NewLen + Order));
  341.     FirstIndex = BspKnotLastIndexL(KnotVector, *Len + Order, t) + 1;
  342.  
  343.     /* Copy all the knot before the knot t. */
  344.     CAGD_GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * FirstIndex);
  345.  
  346.     /* Insert Mult knot of value t. */
  347.     for (i = FirstIndex; i < FirstIndex + Mult; i++)
  348.     NewKnotVector[i] = t;
  349.  
  350.     /* And copy the second part. */
  351.     CAGD_GEN_COPY(&NewKnotVector[FirstIndex + Mult],
  352.           &KnotVector[FirstIndex + CurrentMult],
  353.           sizeof(CagdRType) * (*Len + Order - FirstIndex - Mult + 1));
  354.  
  355.     *Len = NewLen;
  356.     return NewKnotVector;
  357. }
  358.  
  359. /******************************************************************************
  360. * Returns the multiplicity of knot t in knot vector KnotVector, zero if none. *
  361. ******************************************************************************/
  362. int BspKnotFindMult(CagdRType *KnotVector, int Order, int Len, CagdRType t)
  363. {
  364.     int LastIndex, FirstIndex;
  365.  
  366.     if (!BspKnotParamInDomain(KnotVector, Len, Order, t))
  367.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  368.  
  369.     FirstIndex = BspKnotLastIndexL(KnotVector, Len, t) + 1;
  370.  
  371.     for (LastIndex = FirstIndex;
  372.      LastIndex < Len && APX_EQ(KnotVector[LastIndex], t);
  373.      LastIndex++);
  374.  
  375.     return LastIndex - FirstIndex;
  376. }
  377.  
  378. /******************************************************************************
  379. * Scans the given knot vector to a potential C1 discontinuity.              *
  380. *   Returns TRUE if found one and set t to its parameter value.              *
  381. *   Assumes knot vector has open end condition.                      *
  382. ******************************************************************************/
  383. CagdBType BspKnotC1Discont(CagdRType *KnotVector, int Order, int Length,
  384.                                 CagdRType *t)
  385. {
  386.     int i, Count;
  387.     CagdRType
  388.     LastT = KnotVector[0] - 1.0;
  389.  
  390.     for (i = Order, Count = 0; i < Length; i++) {
  391.     if (APX_EQ(LastT, KnotVector[i]))
  392.         Count++;
  393.     else {
  394.         Count = 1;
  395.         LastT = KnotVector[i];
  396.     }
  397.  
  398.     if (Count >= Order - 1) {
  399.         *t = LastT;
  400.         return TRUE;
  401.     }
  402.     }
  403.  
  404.     return FALSE;
  405. }
  406.  
  407. /******************************************************************************
  408. * Scans the given knot vector to aall potential C1 discontinuity. Returns     *
  409. * a vector holding the parameter values of C1 discontinuities, NULL of non    *
  410. * found.                                      *
  411. *   Sets n to length of returned vector.                      *
  412. *   Assumes knot vector has open end condition.                      *
  413. ******************************************************************************/
  414. CagdRType *BspKnotAllC1Discont(CagdRType *KnotVector, int Order, int Length,
  415.                                     int *n)
  416. {
  417.     int i, Count,
  418.     C1DiscontCount = 0;
  419.     CagdRType
  420.     LastT = KnotVector[0] - 1.0,
  421.     *C1Disconts = (CagdRType *) IritMalloc(sizeof(CagdRType) * Length);
  422.  
  423.     for (i = Order, Count = 0; i < Length; i++) {
  424.     if (APX_EQ(LastT, KnotVector[i]))
  425.         Count++;
  426.     else {
  427.         Count = 1;
  428.         LastT = KnotVector[i];
  429.     }
  430.  
  431.     if (Count >= Order - 1) {
  432.         C1Disconts[C1DiscontCount++] = LastT;
  433.         Count = 0;
  434.     }
  435.     }
  436.  
  437.     if ((*n = C1DiscontCount) > 0) {
  438.     return C1Disconts;
  439.     }
  440.     else {
  441.     IritFree((VoidPtr) C1Disconts);
  442.     return NULL;
  443.     }
  444. }
  445.  
  446. /*****************************************************************************
  447. * Routine to determine where to sample along the provided paramtric domain,  *
  448. * given the C1 discontinuities along it.                     *
  449. *   Returns a vector of length NumSamples.                     *
  450. *   If C1Disconts != NULL (NumC1Disconts > 0) this vector is being freed.    *
  451. *****************************************************************************/
  452. CagdRType *BspKnotParamValues(CagdRType PMin, CagdRType PMax, int NumSamples,
  453.                   CagdRType *C1Disconts, int NumC1Disconts)
  454. {
  455.     int i, CrntIndex, NextIndex, CrntDiscont;
  456.     CagdRType Step,
  457.     *Samples = (CagdRType *) IritMalloc(sizeof(CagdRType) * NumSamples);
  458.  
  459.     Samples[0] = PMin;
  460.     if (NumSamples <= 1)
  461.     return Samples;
  462.     Samples[NumSamples - 1] = PMax;
  463.     if (NumSamples <= 2)
  464.     return Samples;
  465.  
  466.     if (NumC1Disconts > NumSamples - 2) {
  467.         /* More C1 discontinuities than we are sampling. Grab a subset of    */
  468.     /* The discontinuities as the sampling set.                 */
  469.         Step = ((CagdRType) (NumC1Disconts - 1)) / (NumSamples - 2) - EPSILON;
  470.     for (i = 0; i < NumSamples - 2; i++)
  471.             Samples[i + 1] = C1Disconts[(int) (i * Step)];
  472.     }
  473.     else {
  474.     /* More samples than C1 discontinuites. Uniformly distribute the C1  */
  475.     /* discontinuites between the samples and linearly interpolate the   */
  476.     /* sample values in between.                         */
  477.         Step = ((CagdRType) (NumC1Disconts + 1)) / (NumSamples - 2);
  478.     CrntIndex = CrntDiscont = 0;
  479.  
  480.     while (CrntDiscont < NumC1Disconts) {
  481.         NextIndex = (int) ((CrntDiscont + 1) / Step);
  482.         Samples[NextIndex] = C1Disconts[CrntDiscont++];
  483.  
  484.         for (i = CrntIndex + 1; i < NextIndex; i++) {
  485.         CagdRType
  486.             t = (i - CrntIndex) / ((CagdRType) (NextIndex - CrntIndex)),
  487.             t1 = 1.0 - t;
  488.  
  489.         Samples[i] = Samples[CrntIndex] * t1 + Samples[NextIndex] * t;
  490.         }
  491.  
  492.         CrntIndex = NextIndex;
  493.     }
  494.  
  495.     /* Interpolate the last interval: */
  496.         for (i = CrntIndex + 1; i < NumSamples - 1; i++) {
  497.             CagdRType
  498.             t = (i - CrntIndex) / ((CagdRType) (NumSamples - 1 - CrntIndex)),
  499.             t1 = 1.0 - t;
  500.  
  501.             Samples[i] = Samples[CrntIndex] * t1 + Samples[NumSamples - 1] * t;
  502.         }
  503.     }
  504.  
  505.     if (C1Disconts != NULL)
  506.     IritFree((VoidPtr) C1Disconts);
  507.  
  508.     return Samples;
  509. }
  510.  
  511. /*****************************************************************************
  512. * Given a knot vector, make sure adjacent knots that are close "enough" are  *
  513. * actually identical. Important for robustness of subdiv/refinement algs.    *
  514. *****************************************************************************/
  515. void BspKnotMakeRobustKV(CagdRType *KV, int Len)
  516. {
  517.     int i;
  518.  
  519.     for (i = 1; i < Len; i++)
  520.     if (KV[i] - KV[i - 1] < KNOT_IS_THE_SAME && KV[i] != KV[i - 1])
  521.         KV[i] = KV[i - 1];
  522. }
  523.